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Abstract 

Plankton is the productive base of aquatic ecosystems and plays a major 
role in the global control of atmospheric carbon dioxide. Nevertheless, after 
intensive study, the factors that drive its spatial distribution are still far from 
being clear. The models proposed so far show very limited agreement with 
actual data as many of their results are not consistent with field observations. 
Here we show that fluctuations and turbulent diffusion in standard prey- 
predator models are able to accurately and consistently explain plankton 
field observations at mesoscales (1-100 km). This includes not only the spatial 
pattern but also its temporal evolution. We explicitly elucidate the interplay 
between physical and biological factors, suggesting that the form in which 
small scale biotic fluctuations are transferred to larger scales may constitute 
one of the key elements in determining the spatial distribution of plankton 
in the sea. 
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Understanding how complex ecosystems work often relies on simplified mod- 
els that disregard many details of the actual system while retaining the es- 
sential information |T], |3fl. In the case of marine ecosystems, not only 
the simplest approaches failed to explain the spatial distribution of plank- 
ton populations but also more sophisticated models were unable to account 
consistently for the most remarkable features fl4|. Even now, it is not clearly 
understood why under apparent homogeneous conditions, such as temper- 
ature and nutrients, plankton is still patchilly distributed. This particular 
situation is of remarkable importance since in the absence of external sources 
of patchiness the pattern must arise as a mere consequence of the interactions 
between the individuals. The typical form of these patterns is illustrated in 
Fig. [I]. The main trait is that zooplankton is more patchilly distributed than 
phytoplankton 

The most intuitive model that can be proposed to explain plankton dy- 
namics [2|] considers the population densities of prey (phytoplankton), N, 
and predators (zooplankton), P: 

^ = F N (N,P) + D N V 2 N , (1) 
BP 

^- = F P (N,P) + D P V 2 P , (2) 

where Dn and Dp are diffusion coefficients; and Fjv and Fp are functions 
that account for the interaction between both species. This class of mod- 
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els is the most frequently used in theories on pattern formation in ecol- 
ogy IS H' 01- I n t ne case °f plankton, they were able to display spatial 
heterogeneity under homogeneous conditions ||. Zooplankton, however, was 
less patchilly distributed than phytoplankton, in contradiction with the ob- 
served pattern |]. 

There are two relevant features that are not taken into account by this 
kind of models. First, diffusion in the sea is not quantitatively well modeled 
by usual Fickian diffusion ||. Both types of diffusion processes will tend 
to spread and mix the populations, but the specific form in which this is 
achieved is different. Second, there is always an intrinsic stochasticity asso- 



ciated with the dynamics of the population (JTy, |TlJ. From birth to death, all 
processes share some degree of chance. The way randomness manifests in the 
dynamics of the individuals depends on the scale we are looking at |TJ; de- 
terministic equations are expected to be valid in the limit of high numbers of 
individuals [p!2| . Therefore, a deterministic description may be a reasonable 
one for phytoplankton alone, but this does not need to be so for zooplankton 
which has much fewer individuals JT3J. More importantly, while phytoplank- 
ton interacts mainly with zooplankton, zooplankton interacts also with fish 
and whales which are far from being evenly distributed. 

These two additional features have been incorporated in a prey-predator 
model: 

dN 

— = F N (N,P)-v-VN, (3) 
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^ = F P (N,P)-v-VP + m, (4) 

where dispersal is given by advection with a velocity field v [= v(r)] that 
depends on the position r, and where a noise term £(t) has been included. 

In general, the effects of the advective terms depend on the precise form 
of the velocity field. For some turbulent fields ||14|| . the effect of advection 
can be simplified as follows: given a passive field f(r,t) which evolves as 



| = (5) 



the spatial Fourier transform of f(r,t) follows from 

-D\kff k , (6) 



__nltlj8- 



dt 

where k is the wave number and D a constant. In this case, advection 
can effectively be viewed as a diffusion process with a diffusion coefficient 
D e ff(k) = D | k | 2-/3 that depends on the scale. In contrast to usual Fickian 
diffusion, the variance of the field is not proportional to t but is given by 
(r 2 ) ~ t 2 ^ . This is the type of time dependence observed for the dispersion 



of tracers in the sea |9], [L5|], from which one can obtain the explicit value of 
the parameter (3. 

To render our model analytically tractable, we consider the system around 
a stable state. Fluctuations in zooplankton, £, move the system away from 
equilibrium. If the fluctuations are not too large, we can perform a linear 
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expansion of F P and F N : 



F N {N,P) = c N -a n N-a 12 P, 

F P (N,P) = c P + a 21 N - a 22 P . (7) 

Here cn, cp, an, a\ 2 , a 2 ±, and 022 are positive constants. For the sim- 



plest form of the noise term ||12|| . Gaussian white and uncorrected in space 
[(£(r, t)) = and (£(f, t)f (7*, t')) = 2o 2 b{f - r)5(t' - t)], and for a 22 ~ the 
variance spectra are given by 

22 2 

S N (k) = —~ a i 2(T - ~ and S P (k) = . (8) 

V ; (D N + Dp)D N D P ^ ! Dp y ' 

where Dn = D\k\@ + an, Dp = D\k\P + a 22 , and a 2 is the intensity of 
the noise source. The assumptions involved do not substantially constrain 
the applicability of the results. When a 2 2 is not negligible, the expressions 
become more involved but the qualitative behavior is still the same. In par- 
ticular, the high wave-number limit remains unchanged. On the other hand, 
the type of noise we have considered is quite general and can arise, among 
others, from a random distribution of predators feeding on zooplankton or 
even from the birth process itself |16|]. Other types of noise with different 



properties, e.g. as those induced by turbulence |14], |4|, are certainly present 
but we assume that they are not relevant for the spectral properties of the 
pattern at the mesoscales. 
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The variance spectra obtained from previous equations display a power- 
law region with exponent —3(3 for the phytoplankton and — (3 for the zoo- 
plankton. The value of (3 ~ 0.87 inferred from diffusion in the sea [To|] 
leads to exponents —2.6 and —0.8, both in excellent agreement with field 
data [[[], |l7], [ll|. It is worth emphasizing that the power law behavior ap- 
pears only for sufficiently high wavenumber (short scales); for low wavenum- 
bers (long scales), the variance spectra is flat, as observed in most field 
data H O, EH- 

In the same way, one can compute the coherence between two patterns at 
different times which provides information about the global dynamics. 



For the phytoplankton this quantity is given by 

N(k, At) D P e 3NAt - D N e 6pAt 



N(k,0) Dp-D 



(9) 



N 



where N(k,At) = J™ \ \N(k,u)\\ 2 e- iujAt du. This result indicates that short 
scales lose their correlation faster than long ones and that eventually the 
whole pattern will be decorrelated, as observed in satellite measurements |18| . 



In Fig. H we plot the typical form of the variance spectra and the squared 
coherence for different time lags. Both of them are in excellent agreement 
with field data |l|, [17], |1| |l9f . Remarkably, the main properties of the pat- 



tern already appear in the linear regime. Therefore, nonlinear interactions 
that drive the system towards a stable state will lead to similar results. To 
study this aspect in more detail, we have performed numerical simulations 
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for typical nonlinear interactions as explained in the caption of Fig. |. The 
resulting two-dimensional spatial distribution, transects, and variance spec- 
tra (shown in Fig. |^) agree with both the linear model and field data. Other 
types of nonlinear interactions — e.g. different functional responses — as 
well as different types of noise — e.g. acting on zooplankton growth rate — 
also produce similar results (data not shown). 

Field observations indicate that the power law region of the variance spec- 
tra and the value of the exponent of this power law are robust properties of 
the system; i.e. these properties are present under a wide variety of con- 
ditions. In our model, there is always a power law region whose exponent 
does not depend on biotic factors but is completely determined by the spe- 
cific form in which turbulent diffusion acts on the system. This provides a 
straightforward explanation of the predominance of the observed exponents 
for the phytoplankton falling between —3 and —2. These are the values that 
arise for 2D {(3 = 1) and 3D (/3 = 2/3) isotropic turbulence, respectively |T^| . 
In the sea, the value of this exponent will depend on the particular situation, 
but it is reasonable to assume that it will be between those of 2D and 3D 
isotropic turbulence, as the available data shows (9|, jljj . There are also non- 
robust properties, such as the region where the variance spectra turns flat. 
In the model, this depends on many factors: e.g. growth and death rate, 
and turbulence. Field data shows that, indeed, the position of this region 
exhibits great variability and that sometimes it is not even present in the 
range of scales observed. 
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Turbulent diffusion and noise are two obvious features that have already 
been considered in the context of marine ecosystems, but none of them by 
itself has been able to explain the mesoscale patterns. In particular, it is 
well known that noise generates variability, i.e. that noise can be a source 



of patchiness |20|, |21|| . For instance, reaction-diffusion prey-predator models 
with noise produce patterns that at a glance strongly resemble those observed 
in the sea (22|, |23)1 . The exponents obtained (—6 for the phytoplankton and 
—2 for the zooplankton), however, are far from the observed ones. This 
quantitative, but not qualitative, disagreement is due to the dependence 
of the effective diffusion coefficient with the scale. Thus, reaction-diffusion 
models are unable to integrate correctly the scale dependence of the physical 
properties of the environment. When this is taken into account, noise not 
only generates patterns but is also able to produce the right ones. 

Turbulence plays a somehow ambivalent role. It can act in the same way 
as diffusion does (transferring variance from smaller to larger scales) and 
also in the opposite way (from larger to smaller scales). These two types of 
processes are referred to as turbulent diffusion and turbulent stirring, respec- 
tively. The former is the one we have considered in our model. It was already 
considered in Ref. fl|| together with the type of noise that turbulence induces 
but without the noise that can arise from biotic factors. The latter only 
plays an important role when some degree of environmental heterogeneity is 



present IQ. Indeed, it has been shown that turbulent stirring can generate 
patterns that resemble the observed ones if spatial heterogeneity and time 
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delays are introduced in the model |25|. The type of time delays introduced, 
however, can lead non-realistic situations such as growing zooplankton in the 
absence of phytoplankton. 

Finally, it is important to emphasize that in our case noise is the key 
element that allows moving from the individual to the population description. 
Our results suggest that zooplankton dynamics at lower scales affects the 
pattern at the mesoscale in the same way as noise does. Considering a more 
detailed description is not necessary to explain and to understand the main 
characteristics of the pattern. This does not mean that the actual dynamics 
of zooplankton is not important at all: its growth rate, its survival, and the 
intensity of noise itself depend, among other factors, on how zooplankton 



aggregates and on how it avoids its predators |26|]. It rather means that 
under a wide range of conditions all those intricate mechanisms will lead to 
patterns with properties as those induced by noise. 
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Figure Captions 

FIGURE 1: Transects obtained from field observations for (a) phytoplank- 
ton and (b) zooplankton (redrawn from Ref. 0). Here Chlorophyll a is a 
measure of phytoplankton activity. 

FIGURE 2: (a) Variance spectra S(k) for prey and predators from Eqs. 
d), (|), and (g). (b) Square coherence [N(k, At)/N(k, 0)] 2 for prey [Eq. @] 
for time lags of 1, 6, and 7 days. The values of the parameters are D = 12, 
an = 0.3, a,22 = 0.05, and [3 = 1. The length and time units are Km and 
days, respectively. 

FIGURE 3: Numerical simulations for a prey-predator model given by 
Eqs. d) and (§, with F N (N, P) = rN (1 - N/K)-cPf(N) and F P (N, P) = 
P{gf{N) - e). The functional response is f(N) = N 2 /(l + N 2 ). Here, K, 
c, g, r and e are positive constants. The velocity field v = v(r) consists 
of a series of vortices distributed as in Ref. [^7]] and the noise is assumed 



to be Gaussian with zero mean and correlation function (f% t)^ (f 7 ', t')) = 
2[aP(r, t)] 2 S(r' — r)5(t' — t). Typical transects for (a) prey and (b) predators, 
(c) Variance spectra for prey and predators. These results were obtained 
by discretizing the corresponding equations on a 250 x 250 two-dimensional 



mesh Pq| , with periodic boundary conditions and then by using a standard 



method for integrating stochastic differential equations p9| . The values of 
the parameters are r = 0.3, K = 4, c = 2, g = 0.1, e = 0.05, and o = 3.5. 
The length and time units are Km and days, respectively. The size of each 
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cell of the discretized mesh is 0.25 x 0.25 Km 2 . 
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